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ABSTRACT 

We have reexamined the similarity solution for a self-gravitating isothermal 
gas sphere and examined implication to star formation in a turbulent cloud. 
When parameters are adequately chosen, the similarity solution expresses an 
accreting isothermal gas sphere bounded by a spherical shock wave. The mass 
and radius of the sphere increases in proportion to the time, while the central 
density decreases in proportion to the inverse square of time. The similarity 
solution is specified by the accretion rate and the infall velocity. The accretion 
rate has an upper limit for a given infall velocity. When the accretion rate is 
below the upper limit, there exist a pair of similarity solutions for a given set of 
the accretion rate and infall velocity. One of them is confirmed to be unstable 
against a spherical perturbation. This means that the gas sphere collapses to 
initiate star formation only when the accretion rate is larger than the upper limit. 
We have also examined stability of the similarity solution against non-spherical 
perturbation. Non-spherical perturbations are found to be damped. 

Subject headings: accretion — hydrodynamics — shock waves — stars: formation 



INTRODUCTION 



Similarity solutions have contributed ver y much t o our unde rstanding of sta r formation 
process. The classical similarity solution by iLarsonl (119691 ) and iPenstonl (119691 ) elucidated 
the runaway nature of gravitational collapse. The density increases in proportion to the 
inverse squa re of t he tim e du ring the runaw ay collapse phase. We learned from the similarity 
solutions of Shul ( 1977 ) and Hunter ( 1977 ) that the accretion rate of a protostar is of the 
order of c s 3 /G where G and c s denote the gravitational constant and the isothermal sound 
speed of gas. 
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The similarity solution is al so used to evaluate the effects of rotati on and magnetic 
field. The similarity solutions of iNarita. Hayashi. fc Miyamal (119841 ) and ISaigo &: Hanawa 
( 119981 ) indicated that the runaway collapse cannot be prevented by rotation if once initi- 
ated. Collapse of a rota t ing m agnetized gas cloud is described by the similarity solution of 



Krasnopolsky fc 



^onigll (120021) . Ambi p olar diffusion is tak en into account in the the simi 



larity solution of lAdams fc Shul (120071). ITsai fc Hsul (119951 ) extended the similarity solution 
to include shock wave. IShu et al.l (120021 ) extended the similarity solution involving a shock 
wave for application to champagne phase of an H II region. 



Tsai &: Hsul (119951 ) found two classes of similarity solution; the first class describes 
accretion onto protostar while the second one does failure of star formation. The central 
density decreases in proportion to the inverse square of the time, p c oc t~ 2 , in the second class 
solution. Although this solution has not gained much attention thus far, it provides an insight 
on dynamical compression of a molecular cloud core. If a dense clump of gas is compressed 
by an external force, the temporal increase in the density may trigger gravitational collapse 
and star formation. One can surmise existence of threshold of gravitational collapse. If the 
dynamical compression is either weak or short, the clump will bounce back to expansion. 
A sho ck wave will be formed when accret ing gas is stoped by the expansion (lAdams fc Shu 
20071 . see, e.g.). The similarity solution of ITsai &: Hsul (119951 ) demonstrated that a spherical 
cloud can expand even when it is steadily compressed by a shock wave. On the other hand, 
the shock compressed gas sphere will collapse owing to its self gravity if the shock is strong 
and lasts for a long enough period. 



In this paper we reexamine the similarity solution of ITsai fc Hsul (119951 ) while keeping 
its negative implication in mind. We find the condition for existence of similarity solution 
describing expansion of a gas sphere. Conversely it will tell us condition for a shock com- 
pressed gas sphere to collapse by its self gravity. We also study stability of the similarity 
solution. The similarity solution denies collapse due to the self gravity only when it is stable. 

We review the similarity solution in §2.1 and show the method of linear stability analysis 
in §2.2. Technical details on the stability analysis are given in Appendix. Properties of 
similarity solutions, such as accretion rate and infall velocity, are shown in §3.1. Stability of 
the similarity solution is given in §3.2 and §3.3. We discuss implications of our analysis in 
$4. 
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2. Model and Methods of Computation 

2.1. Similarity Solution 

We consider an isothermal gas of which distribution is spherically symmetric. Then the 
hydrodynamical equations are expressed as 

dv + v dV + 1 9P + ° Mr (2) 
dt dr p dr r 2 

dMr -A*r*p (3) 



and 



dr 
where 

P = c\p. (4) 

Here, the symbols, p, v, P, M r , G, and c s denote the density, velocity, pressure, mass 
inside the radius r, gravitational constant, and the isothermal sound speed, respectively As 
originally shown by Larson (1969) and Penston (1969), the hydrodynamical equations have 
similarity solutions, 



v{r,t) = H®, (6) 

C s 

M r (r,t) = ^MO (7) 



where 



r 



We restrict ourselves in the case of t > in the following. Substituting Equations (JSj) 
through dSJ) into Equations flTJ through (|2D we obtain 

dQ [/i - 2 (£ - u) 2 } q ^ 



[(« - o 2 - 1] £ 
«9« (/x - 2) (£ - u) 



(10) 



5e [(« - e) 2 - i] e 

and 
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We assume that the density is finite at the center since we are interested in application 
to a starless core, i.e., case of no star formation. Then the density and velocity should be 
expressed as 



Q 



u 





Po 






6 ' 










3 


45 



^ 45 



Qo 



(12) 
(13) 



near the center. Following iTsai fc Hsul (119951 ) we assume that the flow has a shock wave at 
£ = £ s h- Then the Rankine-Hugonio relation gives us the condition, 



Q- 



sh 



U. 



and 



£sh 

(U+ - £sh) (U+ - £sh) 



1 . 



(14) 
(15) 



where g + and £>_ denote the the densities of the pre- and post-shocked gases, respectively, 
andu + and w_ do the velocities of them, respectively. 

In the region far from the origin, a solution of Equations (Q and (flOl) approaches the 
asymptotic solution, 



M 



Q 



u 



^inf£ 2 

- Vin{ + o(r x ) 



(16) 
(17) 



The symbol, v- m t, denotes the infall velocity at the infinity while M denotes the accretion 
rate. The similarity solution can be specified either by (g c , £ s h) or by (tw, M). 

Note that Larson-Penston solution has the same asymptotic form. Only when v- m f and 
M vanish, the solutio n Equations and fllOp have a different asymptotic form, "plus 
solutions" of H Jl977h . 

We integrate equations (jHJ) and (TT0|) by the 4th order Runge-Kutta method from £ = 
for a given set of g c and £ s h to obtain a similarity solution. The infall velocity and accretion 
rate are obtained numerically as a function of q c and £ s h- We also obtain the mass enclosed 
in the shock front, 



M r 



'11 
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2.2. Spherical and Non-spherical Perturbations 

We have performed a normal mode analysis to examine the stability of the similarity 
solution. In the analysis, the density is assumed to be expressed as 

,(r, t) = ' iiitr y 8 - a , ( 19 ) 

where Y™ (9, <p) denotes the spherical harmonic function. This particular form is chosen be- 
cause the similarity solution has no specific timescale. An eigenmode has a growth timescale 
and an unstable perturbation grows exponentially, when the unperturbed state is stationary 
and has a specific timescale. When the physical quantities vary according to a power law 
in time in an un perturbed state, the eig e nmod e grows (or decay) in proportion to a power 



of time, \t\ a . See lHanawa fc Matsumotol (119991 ) for the justification of Equation (fT9l) . They 



analyzed linear stability of the Larson-Penston solution against a non-spherical perturba- 
tion, using the coordinates, (£, 9, <p, In t) instead of the ordinary spherical coordinates, 
(r, 9, if, t). It is shown that the similarity solution can be expressed as a steady state and 
an unstable mode grows in propotion to exp (er In t) = t a in the coodinates. 

The power index, a, is obtained as an eigenvalue of the perturbation equations as shown 
in Appendix. When the real part of a is positive, the mode grows in time and the similarity 
solution is unstable. We call the real part of a, the growth index in the following according 
to the suggestion by F. H. Shu, the referee of this paper. When a is complex, the mode 
grows (or decays) while oscillating. The imaginary part of a is related with the frequency of 
oscillation, although the physical oscillation period increases with time. 



3. Results 

3.1. Similarity Solutions 

First we obtained a series of similarity solutions for a given q c by increasing £ s h- The 
infall velocity, tw, decreases monotonically with increase in £ sh . It reaches t> in f = at some 
£ s h and the similarity solution terminates. By compiling the similarity solutions, we obtained 
M as a function of t>i n f and £ s h- The curves denote M as a function of £ s h for given t> in f in 
Figure [H 

The accretion rate, M, is maximum at a certain £ s h for a given t> in f. It increases with 
increase in £ s h at a small £ s h while it decreases with increase in £ s h at a large £ s h- The former 
is denoted by thin curves while the latter is by thick ones. This means that there exists two 
similarity solutions for a given set of tw and M. 



- 6- 



Figure [2] shows two similarity solutions having f; n f = 1.4 and M = 1.204. The solid 
curves denote the solution of £ s h = 0.986 while the dashed curves do that of £ s h = 0.574. 
Both the solutions have the same density and velocity distributions in the region of £ > 
0.986. The main difference is the location of the shock front. The shock compressed gas 
sphere is denser and expands more slowly in the sol ution denoted by the dashed curve. Since 
the expansion is slow, the solution is similar to the iBonnorl (119561 ) - lEbertl (119551 ) solution for 
a self-gravitating isothermal gas sphere. As well as the Bonnor-Ebert sphere, the latter 
solution is shown to be unstable (§3). 

Figure [3] indicates that the similarity solution has an upper limit on M for a given i^f. 
The upper limit is highest M max = 1-312 at tw = 1.64 . Similarity solutions do not exist 
for M > 1.312. The implication of non-existence is discussed in §4. 



3.2. Spherical Perturbations 

Figure H] shows the growth index of the spherical perturbation, function of £ s h 

for a series of similarity solutions having a given v- m f. The solid curves denote the growth 
index of modes having real index. The dashed lines denote the real part of complex indecies. 
One of the index is positive (a r > 0) and the similarity solution is unstable only when the 
shock radius (£ S h) is smaller than a critical value. The condition of neutral stability coincides 
with that of maximum accretion rate for a given v- m { as expected. When £ s h is a little larger 
than the critical, the solution is stable and the most slowly damping mode has a real index. 
When £ s h is large enough, the solution is stable and the most slowly damping mode has a 
complex index. 

Our survey is limited to modes having low indecies, i.e., in the range < 1.0. It 
is however unlikely that we have missed unstable spherical perturbations. A spherical per- 
turbation induces only sound waves and they are confined within the gas sphere since it is 
bounded by the shock wave. A high frequency sound wave has a shorter wavelength and is 
unlikely to be Jeans unstable. 



3.3. t = 2 Mode 

We have studied 1 = 2 mode as a typical non-spherical perturbation. This is in part 
because the dipole [l = 1) mode is unlikely to be excited. If the dipole mode grows, the 
inner and outer parts of the gas sphere should move in the opposite direction each other to 
keep the center of gravity. 
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Figure [6] denotes the eigenfrequencies of the i = 2 modes for similarity solutions of 
f in f = 1.0. The abscissa denotes £ s h while the ordinate denotes the index, a. All the modes 
are damping (c^ < —1). The solid lines denote the mode having the smallest damping index. 
The mode has a small imaginary part (<7j ~ ±0.1). The mode having the second smallest 
damping index is denoted by the dashed line in Figure [6] and has a real eigenfrequency (pure 
damping). The mode having the third smallest damping index has an imaginary part in 
the eigenfrequency. The imaginary part is similar to that of the mode having the smallest 
damping index. These three modes have a similar damping index of —1.1 < a r < —1.0. 
The other mode (dash-dotted line in Fig. [6]) has a much larger damping index. 

Figure [7| is the same as Figure [6] but for v in f = 4.0. Again, all the modes are damping. 
The damping index is larger than unity (o> < —1). The oscillation frequency of the smallest 
damping mode is larger than those of Vi a t = 1.0. 

One might ask the reason why the similarity solution is stable against a bar (£ — 2) 
mode. We think that the bar mode is damped by expansion of the gas sphere. Since the 
radius of the shock front increases with the time, the asphericity of the shock compressed gas 
sphere decreases unless the displacement grows fa ster than the radius . Whe n a gas sphere is 



collapse, the bar mode can be excite d as sh own by lLin. Mestel. fc Shul (119651 ) for the pressure 



less gas and lHanawa &: Matsumotol (119991 ) for an isothermal gas. 



We have not yet studied the modes of £ > 3. However they are also unlikely to be 
unstable, since the self-gravity does not excite a short wave perturbation. 



4. Discussion 

We obtained the critical accretion rate above which there exits no similarity solution. 
The critical rate can be interpreted as the minimum accretion rate for a high density clump 
to initiate self-gravitational collapse. The critical accretion rate can be rewritten as 



dM 



dt 

for v > 3c s in the dimensional form. 



Equation ([20]) gives us an estimate for a converging flow to initiate gravitational collapse. 
We shall consider a spherical region of which surface is surrounded by a converging flow. The 
radius, inflow velocity, and density are assumed to be r, v and p, respectively. Then the 
gravitational collapse will be initiated when the mass accretion rate exceeds the critical, 

4nr 2 pv > (21) 
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Equation (jzTj) can be rewritten as 

2r 0.3 c s 

V > (22) 

Aj v 

where 

Aj = - 8 . (23) 
V47rGp 

Since Aj denotes the Jeans length, Equation ([22]) means that the effective Jeans length 
reduces in proportion to the inverse of the Mach number. 

The Jeans mass is proportional to the cube of the Jeans length for a given density. Thus 
the effective Jeans mass should reduce to 

M Jieff = Mj l^jf-X (24) 

. -3 

V 



= M] UjtiJ ■ (25) 

This implies that compression of sub Jeans mass clump may result in gravitational collapse 
in the region of flow convergence. Note that the effective Jeans mass is several order of 
magnitude smaller than the classical one when v > 3c s . 

The compression should continue for a certain timescale for a dynamically compressed 
clump to collapse by its self gravity. If we evaluate the minimum timescale to be the effective 
Jeans length divided by the flow velocity, it is shorter than the free-fall timescale by a factor 
of the Mach number squared, 

Aj, c fT f v \ 2 foa\ 

r comp ^ — ^ r s (- . (26) 



The timescale can be translated into the wavelength of perturbation. A compressed 
clump can collapse by the self gravity if the wavelength of velocity perturbation is longer than 
the effective Jeans length. If turbulence contains velocity perturbations of long wavelengths, 
gravitational collapse due to dynamical compression will take place somewhere in the cloud. 
In such case we can expect a number of clumps of which masses are much smaller than the 
classical Jeans mass. 
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We have added comparison of iTsai fc Hsul (119951 ) solution with Bonnor-Ebert solution in 
the revised manuscript according to the comments given to the original manuscript. This 
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A. Linear Stability Analysis 

Since the density is given by Equation f|T9|) in our linear stability analysis, the gravita- 
tional potential should be expressed as 

$(r, t) = c\ hMO + r fy{f-)YT(0, <P)) • (Al) 
The velocity, v = (v r , Vg, v^), is assumed to be expressed as 

v r (r, t) = c s [u(0 + t°5u r Y™(6 } ip)} , (A2) 

v e (r,t) = -^f^^W (A3) 

in the spherical coordiantes. We assume that the perturbation is vanishingly small in the 
region very far from the center. In other words, we restrict ourselves to search for an 
instability due to the internal structure of the shock compressed gas sphere. Thus the flow 
is assumed to have no vorticity, 

(£ + 1) 5u r + {£5u 9 ) = 0. (A5) 

The assumption of no vorticity is based on the Kelvin's circulation theorem. It says the 
circulation, 

T = j> u du (A6) 



is conserved for an isothermal (barotropic) gas (see, e.g., Chapter 6 of IShulll992l ). Since the 



Lagrangian loop C expands, the circulation velocity should decrease w ith the time to conserve 



r, alt hough it can grow with the time if measured at a given £. See lHanawa fc Matsumoto 



( 120001 ) for more details on the mathematical treatment of perturbations having vorticity. 
After some manipulation, the perturbation equations are written as 

{* + 1)1, + ±£&Sf) + = (A7) 
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(* + 2)6f + ^{e[2w5f+(l-w 2 )}} 



^ + |)^ + ^7 + ^ = 0, (A8) 

4w = *r, (A9) 



S 7 = Sg - ^ - 1 J 2 ' 5i> (A10) 



where 



w = u-Z, (All) 

SU <> = <^\c ( W5u r + — + 5 A > ( Al2 ) 

(o- + 1)£ V £o / 

<5/ = Qo8u r + w5g, (A13) 

<5 7 = ^<fy. (A14) 

The perturbation equations do not contain the azimuthal wavenumber, m, and accordingly 
thus the growth rate does not depend on m. This is because the unperturbed state (similarity 
solution) is spherically symmetric. 

We must consider the jump condition at the shock front as well as the boundary condi- 
tions to solve the perturbation equation. The density perturbation is discontinuous between 
the pre- and post-shocked flows and contains change due to shift of the shock front. The lat- 
ter introduces the Dirac's delta function proportional to the shift and the density difference 
between the pre- and post flows. After some manipulation we obtain the jump condition, 



(w+ - wJ) [Sf] - (g + - 



w5u r + — 
go. 



0, (A15) 



(a + 1) [5 7 ] + [Sf] = (A16) 

where the bracket denotes the difference between £ = £ s h ± e. The symbols with the 
subscripts, + and — , denote the values at £ = £ s h ± s, respectively. The perturbation in the 
potential, Sip, is continuous even at £ = £ s h- These jump conditions enable us to connect the 
perturbation in the post-shocked flow (£ < £ s h) and that in the pre-shocked flow (£ > £ s h). 



When the perturbation is spherically symmetric (i = 0), Equations (I A 7ft and flAlOp are 
linearly dependent and we obtain 

(7+1 
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for 



0. Thus we need two independent boundary conditions for 



and three for 



+ 0. 



The spherical perturbation should vanish in the pre-shocked since the unperturbed flow 
is supersonic. Otherwise the perturbation diverges when the mode is unstable. Thus the 
effective outer boundary is set at £ = £ s h and the jump condition is applied there. The other 
boundary condition is set so that the velocity perturbation vanishes at the origin, £ = 
and is proportional to £ near the origin. We have integrated the perturbation equation 
numerically with the Runge-Kutta method from the origin and searched the eigenvalue, a, 
by try and error. 

Non-spherical perturbations should also be regular at the origin and vanishingly small at 
infinity. We calculated asymptotic solutions around the origin and t hose around the infinity 
to exa mine the condition for perturbations to be regular according to lHanawa fc Matsumoto 
( 119991 ) . Some asymptotic solutions are divergent since the perturbation equations are sin- 
gular at the origin and infinity. When £ ^ 0, the perturbation should be expressed by the 
asymptotic solution, 



where 



C 



5g 

dip 



+ 6 



-Ait 1 ' 1 - 
A ( a + 1 



C (£ + 2) 



B 



3 ft -3 |i + 



1 



a 



B 



(A18) 
(A19) 

(A20) 
(A21) 



Here the symbols, A and B, denote arbitrary constants. Equations flA18|) and (1A2CJ mean 
that the density and potential perturbations can be expanded by a series of polynomials 
starting from the £-th power in the Cartesian coordinates. Otherwise the perturbation 
diverges at the origin. The velocity perturbation is also expanded by another series of 
polynomials starting from the (i — l)-th order since it should be balanced with gradient of 
the potential perturbation. 



The other boundary condition is given by the asymptotic solution, 

2D (i + 1) g 



8g 
Su r 
dip 



{a + i + 2)(<t + i + 3)^ + 3 
D(£ + l) 



{a + i + 2)^ + 2 
D 



+i ' 



(A22) 
(A23) 
(A24) 
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for £ ^> 1. Here, the symbol I? denotes an arbitrary constant. Equation (1A24I) means 
that the potential perturbation in the region £ ^> 1 is dominated by the aspherical density 
distribution near the center. In other words the density perturbation is too small to affect 
the gravitational potential. Equation (IA23I) means that the velocity perturbation is induced 
by the gravitational perturbation. Equation (1A22I) means that the density perturbation 
is induced by the velocity perturbation. Thus our boundary conditions gurantee that the 
perturbation is induced by internal change in the flow. 

We integrated the perturbation equations both from the origin (£ = 1CT 2 ) and a very 
large £ (~ 100) with the Runge-Kutta method. We surveyed the eigenvalue, a, by examining 
whether the numerically solutions from both ends satisfy the jump condition at the shock 
front. 

REFERENCES 

Adams, F. C, & Shu, F. H. 2007, ApJ, 671, 497 
Bonnor, W. B. 1956, MNRAS, 116, 351 
Ebert. R. 1955, ZAp, 37, 217 

Hanawa, T., & Matsumoto, T. 1999, ApJ, 521, 703 
Hanawa, T., & Matsumoto, T. 2000, ApJ, 540, 962 
Hunter, C. 1977, ApJ, 218, 834 
Krasnopolsky, R., & Kdnigl, A. 2002, ApJ, 580, 987 
Larson, R. B. 1969, MNRAS, 145, 271 

Larson, R. B. 2003, Reports on Progress in Physics, 66, 1651 

Lin, C. C, Mestel, L., & Shu, F. H. 1965, ApJ, 143, 1431 

Narita, S., Hayashi, C, Miyama, S. M. 1984, Prog. Theor. Phys., 72, 1118 

Penston, M. V. 1969, MNRAS, 144, 425 

Saigo, K., & Hanawa, T. 1998, ApJ, 493, 342 

Shu, F. H. 1977, ApJ, 214, 488 



-13- 

Shu, F. H. 1992, The Physics of Astrophysics Volume II Gas Dynamics, (Mill Valley, Uni- 
versity Science Books) 

Shu, F. H., Lizano, S., Galli, D., Canto, J. & Laughlin, G. 2002, ApJ, 580, 969 

Tsai, J. C, Hsu, J. J. L. ApJ, 448, 774 



This preprint was prepared with the AAS IATgX macros v5.2. 




Fig. 1. — The accretion rate, M, is shown as a function of xi^ for a given t> inf . 
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Fig. 2. — The density and velocity distributions are shown for two similarity solutions having 
the same f; n f and dM/dt. 
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Fig. 3. — The critical accretion rate is shown as a function of v- mi . 



-17- 




Fig. 4. — Each curve denotes the real part of eigenfrequency, ay, of a spherical perturbation 
as a function of the shock radius, £ s h for a given f inf . 
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Fig. 5. — The same as Fig. H]but for the imaginary part, cij. 
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Fig. 6. — The eigenfrequency, a, is shown as a function of £ s h for similarity solutions having 
fi n f = 1.0. The solid curves denote the real part of a for non-spherical perturbations having 
i = 2, while dashed curves do that of imaginary part. Mode a has complex eigenfrequencies 
while modes b and c have real eigenfrequencies. 
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Fig. 7. — The same as Fig, [6] but for t> in f = 4. 



